## 

rm(list = ls())

library(tidyverse)
library(pbapply)
library(rdrobust)
library(fixest)
library(broom)

## Def covars

covs <- c('turnout_party_2009_btw', 'soz_vers_beschaeftigte_share',
          'pop_density_km2', 'migration_out_share')

## Get federal election date 

bt <- readRDS("data/data_federal.rds") %>%
  filter(!is.na(treated)) 

## List of outcomes

outcomes <- c("turnout_party", 
              "agg_left_party", 
              "agg_center_party",
              "agg_right_party")

## Source helper function to tidy rdrobust output

source("code/tidy_rd.R")

## Drop the 09 obs

bt <- bt %>% 
  filter(year > 2012)

## First differences

diff_df <- pblapply(outcomes, function(o) {
  out <- bt %>%
    filter(year > 2012) %>%
    pivot_wider(values_from = o, names_from = 'year', id_cols = 'ags',
                names_prefix = 'o') %>%
    mutate(diff = o2017  - o2013) %>%  dplyr::select(ags, diff) 
  ## Rename
  colnames(out)[2] <- o
  
  ## Return this
  out
}) %>%
  reduce(left_join) %>%
  left_join(bt %>% dplyr::select(ags, pop_dec_09, applies_census,
                                 one_of(covs), state_id) %>%
              distinct(ags, .keep_all = T)) %>%
  mutate(runvar = (pop_dec_09 * -1) + 10000) %>%
  filter(applies_census == 1) %>% 
  mutate(east = ifelse(as.numeric(state_id) > 11, 1, 0))

## Make separate East / West DFs

list_dfs <- list(diff_df %>% filter(east == 1), diff_df %>% filter(!east == 1))
list_labels <- c('East', 'West')

## Estimate

out_rd_opt <- pblapply(outcomes, function(o) {
  lapply(1:length(list_dfs), function(i) {
    
    df_use <- list_dfs[[i]]
    df_label <- list_labels[i]
    
    ## Estimate
    
    out <- rdrobust(y = df_use %>% pull(!!o), 
                    x = df_use$runvar, c = 0,
                    covs = df_use[, covs])
    
    ## Return this : Robust B-C SE
    out2 <- out %>% 
      tidy_rd(se_nr = 3) %>% 
      mutate(outcome = o, method = 'robust bias-corrected') %>% 
      dplyr::rename(se = std.error, pval = p.value, 
                    bw = bw_left_h,
                    bw_bias = bw_left_b) %>% 
      mutate(conf.low90 = estimate - qnorm(0.95)*se,
             conf.high90 = estimate + qnorm(0.95)*se,
             sample = df_label)
    
    ## Both
    
    out2
  }) %>% reduce(rbind)
}) %>% reduce(rbind)

## Check

out_rd_opt %>%
  mutate(estimate = round(estimate, 2),
         pval = round(pval, 2)) %>% 
  dplyr::select(outcome, estimate, pval, sample)

##

out_rd_opt <- out_rd_opt %>%
  mutate(outcome = dplyr::recode(outcome, 
                                 `agg_left_party` = 'Left-wing\nparties',
                                 `agg_center_party` = 'Centrist\nparties',
                                 `agg_right_party` = 'Right-wing\nparties',
                                 `turnout_party` = 'Turnout')) %>% 
  mutate(sample = dplyr::recode(sample, 
                                `East` = 'East Germany',
                                `West` = 'West Germany')) %>% 
  mutate(outcome = factor(outcome, levels = unique(outcome)[c(2:4, 1)]))

# Figure A.12: East and West Germany ----

pd <- position_dodge(0.4)

p1 <- ggplot(out_rd_opt %>%
               filter(method != 'regular'), 
             aes(outcome, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted',
             color = 'grey40')  +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.65) +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90), 
                width = 0, size = 1.05) +
  geom_point(shape = 21, fill = 'white', size = 2,
             position = pd) +
  xlab('') + ylab('RD estimate (percentage points)') + theme_bw() +
  facet_wrap(~sample) +
  theme(legend.position = 'bottom')
p1
